#find x to maximized f
#f = (1+bx)^p(1-x)^q
#ifelse(g,bx,-x),P(g==T) = p, P(g==F) = q
n = 100
p = 0.7 #
b = 0.8 #return from investing 1 unit, odds b to 1
kf = (b*p-(1-p))/b #Kelly formula, the fraction of total capital
if(kf < 0){
kf = 0
warning("warning: not to invest ")
}
x = 0.02
er = 1+x*(p*b+p-1)
#(1+b*x)*p+(1-x)*(1-p)
#p+b*p*x+1-x-p+x*p = 1-x+p*x*b+p*x
vr2 = x^2*(1+b)^2*p*(1-p)
#(1+b*x - 1+x-p*x*b-p*x)^2 * p +(1-x - 1+x-p*x*b-p*x)^2*(1-p)
#((1-p)*x*(b+1))^2*p + (p*x*(1+b))^2*(1-p)
sdr = x*(1+b)*sqrt(p*(1-p)) # linear to x
assumed you invest 20% of your capital. then after one period, your expected capital is 1.052 with standard deviation 0.1649727
#the function kelly try to optimized -- the accumulated return.
# as the number of bets increase, W/N -> p, and L/N -> 1-p. so it converge to this function
zf = function(r,b=1.6,p=0.6){(1+b*r)^p*(1-r)^(1-p)}
#in term of the tree method, for n bets
zf.tree = function(x,b=1.6,p=0.6,n=100){
ss = c(1+b*x,1-x)
ps = dbinom(0:n,n,p)
A = c(1,rep(ss[1],n)) %>% cumprod()
B = c(1,rep(ss[2],n)) %>% cumprod() %>% rev()
return(list(ps = ps, values = A*B))
}
#in term of the tree method, for n bets
zft = function(x,b=1.6,p=0.6,n=100){
ans = zf.tree(x,b,p,n)
ps = ans$ps
Vs = ans$values
e.v = sum(ps*Vs)
e.s2 = sum(ps*(Vs-e.v)^2)
return(list(e.ret = e.v, e.sd = sqrt(e.s2)))
}
#in term of the tree method, for n bets
zft.log = function(x,b=1.6,p=0.6,n=100){
ans = zf.tree(x,b,p,n)
ps = ans$ps
Vs = log(ans$values) #use log transform, because the upper end is too big. take it as utility function
e.v = sum(ps*Vs)
e.s2 = sum(ps*(Vs-e.v)^2)
return(list(e.ret = e.v, e.sd = sqrt(e.s2)))
}
zft.mr = function(x,b=1.6,p=0.6,n=100,a=0.95){
ans = zf.tree(x,b,p,n)
psc = cumsum(ans$ps)
e.v = ans$values[which.max(ans$ps)]
a1 = (1-a)/2
a2 = 1-a1
ind1 = which.min(abs(psc -a1))
ind2 = which.min(abs(psc -a2))
e.sd = ans$value[ind2]-ans$value[ind1]
return(list(e.ret = e.v, e.sd = e.sd, e.sd_low = ans$value[ind1], e.sd_hi = ans$value[ind2]))
}
zft.mm = function(x,b=1.6,p=0.6,n=100,w=0.5){
ans = zf.tree(x,b,p,n)
ind = which.max(ans$ps)
e.v = ans$values[ind]
c1 = 0
ind1 = ind-1
while((ind1>0) & (c1 < w/2)){
c1 = c1+ans$ps[ind1]
ind1 = ind1-1
}
c2 = 0
ind2 = ind+1
while((ind2<n+1) & (c2 < w/2)){
c2 = c2+ans$ps[ind2]
ind2 = ind2+1
}
if(ind1 == 0) ind1 = ind1+1
if(ind2 >n) ind2 = ind2-1
e.sd = ans$value[ind2]-ans$value[ind1]
#print(c(ind,ans$ps[ind],ind1,ind2,ans$ps[ind1:ind2]))
return(list(e.ret = e.v, e.sd = e.sd, e.sd_low = ans$value[ind1], e.sd_hi = ans$value[ind2]))
}
#x is the fraction to be used for the bet
gdx_data = function(p,b,x,m_paths=500,n_steps = 100){
ps = c(p,1-p)
ss = c(1+b*x,1-x)
r.ind = sample(c(1,2),m_paths*n_steps,replace = T, prob = ps)
r = ss[r.ind]
gd = data.frame(path = gl(m_paths,n_steps),d.ret = r,o = rep(1:n_steps,m_paths))
gd = gd %>% group_by(path) %>% mutate(c.ret = cumprod(d.ret))
gdx = rbind(data.frame(path = gl(m_paths,1),
d.ret = rep(1,m_paths),
o = rep(0,m_paths),
c.ret = rep(1,m_paths)),gd)
return(gdx)
}
data = gdx_data(p,b,kf,m_paths = 800)
data$x = kf
xs = c(0.02,0.2,0.5,0.7,0.8)
for(x in xs){
gdx = gdx_data(p,b,x)
gdx$x = x
data = rbind(data,gdx)
}
gp = ggplot(data = data,aes(x = o, y = c.ret, group = path)) +
geom_line(alpha = 0.02) +scale_y_log10()+facet_wrap(~x)
print(gp)
x = 0.2
gf = function(x){zf.tree(x,b,p) %>% as.data.frame() %>% mutate(x = x, pc = cumsum(ps),cf95 = pc > 0.025 & pc <0.975 )}
g = gf(x)
ps = g$ps
ind = which.min(abs(cumsum(ps)-0.975))
p.level = ps[ind]
xs = c(0.02,0.2,0.5,0.7,0.8,kf)
g.data = map_dfr(xs,gf)
pg = ggplot(g.data %>% filter(values > 1e-15 & values < 1e6),aes(x = values, y = ps,color = factor(x)))+
geom_point(alpha = 0.7,size = 0.5)+geom_line(aes(group = factor(x)),alpha = 0.5)+
geom_vline(xintercept = 1,linetype="dashed")+
geom_hline(yintercept = p.level,color="red",alpha = 0.4,linetype="dashed")+
#facet_wrap(~factor(x),ncol = 1)+
scale_x_log10()+theme(legend.position = "none")
print(pg)
data.end = data %>% filter(o == max(o)) %>% mutate(log.c.ret = log(c.ret))
data.end.stat = data.end %>% group_by(x) %>%
summarise(e.ret = mean(c.ret),s.ret = sd(c.ret),
m.ret = quantile(c.ret,0.5),
q05.ret = quantile(c.ret,0.05),
q95.ret = quantile(c.ret,0.95),
e.elog.ret = mean(log.c.ret)%>%exp(),s.elog.ret = sd(log.c.ret)%>%exp(),
.groups = "keep")
knitr::kable(data.end.stat,digits = 3)
| x | e.ret | s.ret | m.ret | q05.ret | q95.ret | e.elog.ret | s.elog.ret |
|---|---|---|---|---|---|---|---|
| 0.020 | 1.679 | 0.278 | 1.657 | 1.242 | 2.137 | 1.656 | 1.180 |
| 0.200 | 129.715 | 288.974 | 40.249 | 2.940 | 786.508 | 36.624 | 5.153 |
| 0.325 | 1491.265 | 6923.664 | 80.362 | 0.545 | 6346.273 | 77.107 | 15.409 |
| 0.500 | 21932.848 | 122078.999 | 15.778 | 0.011 | 59611.065 | 21.489 | 121.295 |
| 0.700 | 2614.954 | 32623.556 | 0.007 | 0.000 | 698.848 | 0.003 | 2268.210 |
| 0.800 | 666.551 | 8393.498 | 0.000 | 0.000 | 2.929 | 0.000 | 20291.095 |
gp = data.end %>%
ggplot(aes(x = x, y = c.ret)) +
stat_summary(alpha = 0.8,geom = "bar",fun = "mean",width = 0.03)+
stat_summary(geom = "errorbar",fun.data = "mean_se",width = 0.01)+
geom_hline(yintercept = 1,linetype = "dashed")+
scale_y_log10()
print(gp)
#p = 0.7 #
#b = 0.7 #return from investing 1 unit, odds b to 1
#er = b*p-1*(1-p) #expected return
#kf = (b*p-(1-p))/b
xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.95,0.01)
zd = data.frame(fraction = zx, expected.return = zf(zx,b,p),sdr = zx*(1+b)*sqrt(p*(1-p))/10)
zd2 = data.frame(fraction = xs, expected.return = zf(xs,b,p),sdr = xs*(1+b)*sqrt(p*(1-p))/10)
zd$color = "grey"
zd2$color = "blue"
zd = rbind(zd,zd2)
pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
geom_line(aes(x = fraction, y = expected.return+sdr),color = "grey")+
geom_line(aes(x = fraction, y = expected.return-sdr),color = "grey")+
geom_point(aes(x = kf,y = zf(kf,b,p)),color = "red")+
scale_colour_identity()
print(pg)
#p = 0.7 #
#b = 0.7 #return from investing 1 unit, odds b to 1
#er = b*p-1*(1-p) #expected return
#kf = (b*p-(1-p))/b
xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.8,0.01)
f.zft = function(x){zft.log(x,b,p) }
zd = map_df(zx,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>%
mutate(fraction = zx,sdr = sdr,color = "grey")
zd2 = map_df(xs,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>%
mutate(fraction = xs,sdr = sdr,color = "blue")
zd = rbind(zd,zd2)
pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
geom_line(aes(x = fraction, y = expected.return+sdr),color = "grey")+
geom_line(aes(x = fraction, y = expected.return-sdr),color = "grey")+
geom_point(aes(x = kf,y = f.zft(kf)$e.ret),color = "red")+
scale_colour_identity()+ylab("expected.log.return")
print(pg)
pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return %>% exp(),color = color))+
geom_line(aes(x = fraction, y = (expected.return+sdr/10)%>% exp()),color = "grey")+
geom_line(aes(x = fraction, y = (expected.return-sdr/10)%>% exp()),color = "grey")+
geom_point(aes(x = kf,y = f.zft(kf)$e.ret %>% exp()),color = "red")+
scale_colour_identity()+ylab("expected.log.return")
print(pg)
#p = 0.7 #
#b = 0.7 #return from investing 1 unit, odds b to 1
#kf = (b*p-(1-p))/b
#if(kf < 0) kf = 0
xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.3)} #the central 30% confident zone
zd = map_df(zx,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>%
mutate(fraction = zx,sdr = log1p(sdr),color = "grey")
zd2 = map_df(xs,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>%
mutate(fraction = xs,sdr = log1p(sdr),color = "blue")
zd = rbind(zd,zd2)
pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
geom_line(aes(x = fraction, y = e.sd_low),color = "grey")+
geom_line(aes(x = fraction, y = e.sd_hi),color = "grey")+
geom_point(aes(x = kf,y = f.zft(kf)$e.ret),color = "red")+
scale_colour_identity()+ylab("most.likely.return")
print(pg+scale_y_log10())
print(pg)
zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft(x,b,p)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
scale_y_log10()+scale_x_log10()
print(pg)
zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.log(x,b,p)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")
print(pg)
zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.95)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")
print(pg)
zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.3)}
d1 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d1$a = 0.3
f.zft = function(x){zft.mr(x,b,p,a = 0.6)}
d2 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d2$a = 0.6
f.zft = function(x){zft.mr(x,b,p,a = 0.1)}
d3 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d3$a = 0.1
d = rbind(d1,d2,d3)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
geom_path(aes(x = e.sd,y = e.ret,linetype = factor(a)),color = "grey")
print(pg)
zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mm(x,b,p,w = 0.9)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")
print(pg)
zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mm(x,b,p,w = 0.3)}
d1 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d1$w = 0.3
f.zft = function(x){zft.mm(x,b,p,w = 0.6)}
d2 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d2$w = 0.6
f.zft = function(x){zft.mm(x,b,p,w = 0.1)}
d3 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d3$w = 0.1
d = rbind(d1,d2,d3)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
geom_path(aes(x = e.sd,y = e.ret,linetype = factor(w)),color = "grey")
print(pg)
DiagrammeR::grViz("digraph {
node [shape = oval]
a [label = '1']
b [label = '1+bx']
c [label = '1-x']
a->b [label='p']
a->c [label='1-p']
d [label = 'x']
e [label = '(1+b)x']
f [label = '0']
d->e [label='p']
d->f [label='1-p']
{rank=same;a;d;}
}",
height = 200)
N = 1e3
#pb.data = data.frame(p = runif(N),b = rnorm(N,mean = 1,sd = 0.3)) %>% filter(b >0)
pb.data = data.frame(p = runif(N),b = runif(N,min = 0,max = 2))%>% filter(b >0)
pb.data = pb.data %>% mutate(kf = (b*p-(1-p))/b,jd = kf > 0,kf_ = ifelse(kf>0,kf,0))
pg = ggplot(pb.data)+geom_point(aes(x = p, y = b,shape = jd,color = kf_))+
geom_line(data = data.frame(x = seq(0.33,1,0.01)),aes(x = x,y = 1/x-1),color = "red")+
guides(shape=guide_legend(title="kf > 0"))+
scale_shape_manual(values=c(1,3))
print(pg)